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1 Introduction 



In this Letter we introduce new classes of energy-preserving numerical meth- 
ods for multidimensional canonical Hamiltonian systems by a modification 
of classical discrete gradient schemes in a locally exact way. Conservative 
integrators have been used for many years for simulating dynamics of sys- 
tems of particles [H [2] . More recently discrete gradient methods have been 
developed in the context of geometric numerical integration [3], see [HE]. In 
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particular, numerical integrators constructed by Quispel and his co-workers 
preserve integrals of motion of a given system of ordinary differential equa- 
tions [61 [7J H]. It is well known that preservation of geometric properties by 
numerical algorithms is of considerable advantage [9]. 

A numerical scheme is called locally exact if its linearization coincides 
with the exact discretization of linearized equations PHI El H2]- We point 
out that exact discretization gives the exact solution at each time step. Thus 
two known procedures are combined: approximation of nonlinear systems 
by linear equations and explicit exact discretizations of linear ordinary dif- 
ferential equations with constant coefficients. Exact discretization of the 
linearized equation (i.e., the exponential Euler method) has been considered 
a long time ago [13], compare [TJ]. Recently, a similar notion (preservation 
of the linearization at all fixed points) appeared in [15J . see also [16]. Locally 
exact numerical schemes can be considered as a special case of exponential 
integrators [UJ HS1 US] . Our approach considers non-standard modifications 
of numerical schemes (see [20]), parameterized by some functions (e.g., a ma- 
trix denoted by S). In contrast to the original Mickens approach, where such 
functions are postulated or guessed, we determine them by requiring local 
exactness. The final formulae are sometimes reminiscent of results generated 
by Gautschi-type methods [2T| |22]. 

Locally exact modifications of discrete gradient methods for Hamiltonian 
systems with one degree of freedom were studied in [TUl HI]. We modified 
the discrete gradient scheme in a locally exact way preserving its main ge- 
ometric property: the exact conservation of the energy integral. Numerical 
experiments show that locally exact modifications can increase the accuracy 
by several orders of magnitude. The multidimensional case is much more 
complicated. First results, confined to the coordinate increment discrete 
gradient and its symmetric modification, were reported in [12] . In this paper 
we present general results that are valid for any discrete gradient. 

2 Local exactness 

We consider an ordinary differential equation (ODE) y = F(y) with solution 
y(t) G W 71 (for some m G N and initial condition y(to) = y ), and a difference 
equation with variable time step h n and solution (y n ) C M m . We denote 
h n = t n+ i — t n , hence we get a sequence t n for n G N U {0}. In other words, 
there is a map that takes the sequence (t n ) to the sequence (y n ) and each 



vector y n is the approximate solution at time t n . We say that the difference 
equation is an exact discretization of the ODE if y n = y(t n ) for n G NU {0}. 
It is well known that any linear ODE with constant coefficients admits the 
exact discretization in an explicit form [23] , see also J20J I2H [22]. Moreover, 
exact discretizations were applied in the numerical solution of the classical 
Kepler problem [2SI EZl EE] and the wave equation [28J. The central topic 
of our paper is another fruitful direction of using exact integrators, namely 
the so called locally exact discretizations [10J [28] . Motivated by the results 
of [TOj E] we propose the following definition (see also [T2]). 

Definition 2.1. A numerical scheme y n+1 = Q(y n ,h n ) for an autonomous 
ordinary differential equation y = F(y) is locally exact if there exist a se- 
quence (y n ) such that y n — y n = 0{h n ) and the linearization of the scheme 
around y n is identical with the exact discretization of the differential equation 
linearized around y n (for any n). 

In particular cases we usually assume y n = y n or y n = \ iy n + Vn+i), 
see [TOl El H2]. It is worthwhile to point out that all proofs included in this 
Letter are valid for any choice of (jj n ). 

Similar concept ("linearization-preserving" schemes) has been recently 
formulated in [15]. A scheme is said to be linearization-preserving if it is 
linearization preserving at all fixed points. All locally exact schemes are 
linearization-preserving but, in general, a linearization-preserving scheme 
may not be locally exact. 

If y is near y~ n , then the equation y = F(y) can be approximated by 

u = F'{y n )u + F{y n ) (2.1) 

where F' is the Frechet derivative (Jacobi matrix) of F and 

v = y-Vn- (2.2) 

The exact discretization of the approximated equation (12.1 p is given by 

i/n+i = e h " F '^u n + h nVl (h n F'(y n ))F(y n ) (2.3) 

(compare [121 [13]), where cpi is an analytic function defined by 

°° z k ~ 1 

k=i 



see, e.g., [29] . For z^Owe have ipi(z) = ie z — l)/z. Moreover, 



fe=i 



(2*)! 



where 5^ are Bernoulli numbers (B\ = |, _B 2 = ^j, -B 3 = ^, etc.), compare 
e.g. [5U]. This series converges for |^| < 2-k. Note that 

ehnF'Wn) = j + h n F'(jj n )<p 1 (h n F J (y n )), (2.5) 

which holds also for singular F'(y n ). In this Letter we do not need to assume 
detF'(y n )^0. 

Here and in the next sections we often use analytic functions of matrices, 
see e.g. [3TJ EZ|- If g '■ C — > C is analytic on the disc {z e C : \z\ < r} for 
some r, then it has a locally defined Taylor series g(z) = 'Y^ =a g n z n and we 
can define g(A) by 

oo 

g(A) = J29nA n , (2.6) 

n=0 

for any matrix A (or a continuous linear operator in a Banach space) such 
that ||v4|| ^ r, where || ■ || is a norm satisfying ||AB|| ^ ||A||||I?|| for any 
operators A and B. If all eigenvalues Ai, . . . , A m of a matrix A are disctinct, 
then A has the spectral factorization A = Tdiag(Ai, . . . , A m )T _1 , where fcth 
column of the matrix T is the eigenvector associated with A&, and we can use 
a compact formula g(A) = Tdiag(g(Ai), . . . , <7(A m ))T _1 . 

Locally exact discretizations have some qualitative advantages which fol- 
low directly from their definition. Namely, they preserve all fixed points of 
the considered differential equation and they are linearly stable and A-stable. 
Indeed, locally exact schemes applied to linear equations produce exact so- 
lutions and exact trajectories. 



3 Discrete gradients 

We consider arbitrary multidimensional Hamiltonian systems in canonical 
coordinates: 

. j. dH . OH . 

x =dpk' P= ~d^> (* = !,...,"»), (3.1) 
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where H = H(x,p) and x := (x 1 , . . . ,x m ), p := (p 1 , . . . ,p m ). Denoting 
y := {x,p) T G M? m we obtain a more concise form of the Hamiltonian system: 

y = F(y) , F(y) = SVH , S = ( [ _ 3 / H . (3.2) 



where 



V-ff" = — — = H y := (H y i,H y 2, ... , H y 2 m ) , Hyj = — — . (3.3) 



It is well known that Hamiltonian if is a first integral of the system (13.21) 



since 

£ = E^<™l*> = <V*l™> = 0, (3.4) 

fc=l W 

where the bracket (denoting a scalar product) is defined by the second equal- 
ity. The last equality of f !3.4p holds for any skew- symmetric matrix S. 

Definition 3.1. A discrete gradient 

- (AH AH AH \ 

* \WAy*''"'Ay»") ^ 5) 

is a continuous function VH : M 2m x R 2m 3 (y n ,Vn+i) h ^ Vi?(j/„,y„ + i) G 
R 2m , such that ]3 \Bj: 

{VH(y n ,y n+l ) | y n+ i - y n ) = H(y n+1 ) - H(y n ) , (3.6) 

limVH(y,y) = VH(y) . (3.7) 

Discrete gradients are non- unique, compare JUEIE]- Many functions sat- 
isfy these conditions, for example the so called coordinate increment discrete 
gradient [I]: 

AH H{y* y 2 n , y 3 n , . . . , y 2 n m ) - H(yl y 2 , y 3 n , . . . , y 2 n m ) 



Ay 1 
AH 


** [Vn+li Vn+l> Vn) • 


Vn+l Vn 

..,y^)-H(y' n+l ,yl.. 


■ , y 2 n m ) 


Ay 2 




Vn+l ~ Vn 




AH 


_ H(yl +1 ,y 2 +1 ,.. 


) Vn+l) ~ " yVn+li Vn+li ' 


■■>y 2 n m ) 



(3.8) 



^i/ 4/n+l i/n 



In fact, any permutation of 2m coordinates x ,p> can be identified with 
components oiy. Thus we obtain (2m)! discrete gradients of this type (some 
of them may happen to be identical). 

Having any discrete gradient we can define the related symmetric discrete 
gradient 

V s H{y n ,y n+1 ) = - (VH(y n ,y n+1 ) + VH(y n+1 ,y n )) . (3.9) 

One can easily verify that V S H satisfies conditions ( 13. 6p . ( 13. 7p provided that 
they are satisfied by VH . 

In order to construct locally exact modifications we need to linearize 
discrete gradients defined by ( 13. 8 p and ( 13. 9p . 

Definition 3.2. We define 2m x 2m matrices A,B as follows 

dVH(y,y) 



A= __ dWH(y,y) 



dy 



B 



y = y n dy 



V = Vr, 



(3.10) 



y = Vn y = Vn 

Thus linearization of the discrete gradient around y n is given by 

VH{y n ,y n+1 ) ^H y + Au n+1 + Bv n , (3.11) 

where we denoted v n =y n —y n and H y is evaluated at y n . 

In the case of the coordinate increment discrete gradient (J3.8P entries of 
matrices A, B can be obtained by straightforward calculation (compare [T2]): 

r H yJyk (j>k) ( o (j > k) 

Ajk = < 2^v k v k C? = ^) ' Bjk = \ 2^y k y k U = ^) • (3.12) 

( (j <k) [ H yJyk (j < k) 

Theorem 3.3. For any discrete gradient we have: 



B = A T , A + B = Hyy , 

A + A = Hyy , B + B = Hyy , 

where Hyy is the Hessian matrix of H , evaluated at y r 



(3.13) 



Proof: We expand (|3.6p in a Taylor series, expanding y n and y n +\ around y n . 
Taking into account (|3.1ip . 

H{y n ) = H+{H y \v n ) + -{v n \H m v n ) + ... , (3.14) 

and the analogical expansion for H(y n +i), we observe that the linear parts of 
both sides of (|3.6p are obviously identical. Considering quadratic terms we get the 
following equations: 

(i/ n+1 1 Au n+1 ) = 2< v n+l I H yy v n+1 ), (u n | Bv n ) = -{v n | H yy u n ), (3.15) 

(i/ n+ i | Bv n ) - (Av n+1 1 v n ) = , (3.16) 

which have to be satisfied for any u n , v n+ \. After elementary algebraic manipula- 
tions we obtain: 

(u n+1 1 (A - \H yy ) v n+l ) = , (u n | (B - \Hyy) v n ) = , (3.17) 



{v n+1 \ (B - A T ) v n ) = . (3.18) 

Equations (|3.17p imply that matrices A — \H yy and B — \H yy are skew- 
symmetric. Hence A + A = H yy = B + B. Finally, from (13. 18f) we get B = A 
and, as a consequence, A + B = -f/yy. □ 

Corollary 3.4. Linearization of any symmetric discrete gradient is given by 

V s H(y n ,y n+i ) ^ H y + -H yy (u n + u n+1 ) . (3.19) 

Indeed, in the symmetric case B = A. Therefore we can simplify (13. lip 
using A + B = H yy . 

4 Modified discrete gradient schemes 

The standard discrete gradient scheme is given by 

Vn+l ~ Vn = KSVH , (4.1) 

where VH = VH(y n ,y n+ i) is a given discrete gradient. Replacing S by a 
skew symmetric matrix, we obtain a large class of energy-preserving numer- 
ical schemes, closely related to the discrete gradient method. 



Lemma 4.1. Let A be 2m x 2m matrix (depending on h n ,y n and y n +i) 
such that 



T • ,. 1 . „ „ ( J 



A J = -A , lim — A = S , S=( " : (4.2) 

and Vi? = VH(y n ,y n+ i) is any discrete gradient. Then, the numerical 
scheme 

Vn+i ~Vn = AVH(y n ,y n+1 ) , (4.3) 

zs a consistent energy-preserving integrator for Hamiltonian system A3.2\) . 

Proof: The consistency follows immediately from (I4.2p . The energy preservation 
can be shown in the standard way. Using the scalar product, we multiply both 
sides of ffi2) by VH 



(Vi? | y n+1 - y n ) = (VH \ AVH) . (4.4) 

By (|3.6p the left-hand side equals H(y n+ i) — H(y n ). The right hand side vanishes 
due to the skew symmetry of A. Hence H(y n+ i) = H(y n ). □ 

The following theorem characterizes locally exact modifications of the 
discrete gradient method applied to equation (I3.2p . Any (y n ) satisfying 
Vn — Vn = 0(h n ) is acceptable. 

Theorem 4.2. Numerical scheme y n +i —Vn = A.VH (where A depends on 
y n and h n ) is locally exact if 

A = h n $ 1 S(I + h n A$iS)-\ (4.5) 

where <3>i = ipi(h n F'), F' = SHyy (evaluated aty n ), A is given by li3.10\) and 
we assume that (J + h n A§ 1 S)~ exists. 

Proof: By virtue of (I3.1fp linearization of (14. 3p (around y n ) is given by 

v n +i ~v n = A(Ai/ n+ i + Bv n ) + AHy . (4.6) 

Hence, taking into account SH y = F (where F = F(y n )), we have 

(/ - A^) u n+1 = (/ + AB) v n + AS~ X F . (4.7) 
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The scheme (|4.3p is locally exact if and only if (|4,7j) coincides with f|2.3|) . Therefore, 
we require that 

(/ - AA) e hnF ' =I + AB , (4.8) 



h n (I - AA) $iF = A5 _1 F . (4.9) 

Using Theorem 13.31 (namely, A + B = S~ 1 F') and equation (J2.5J) . we transform 
equation fj4.8|) into 

/i n (J - AA) $iF' = AS^F' . (4.10) 
Both equations (|4.9|) and (|4.1U|) are simultaneously satisifed if 

K (I - AA) $i = AS' 1 (4.11) 
(for invertible F' this sufficient condition is also necessary). Therefore, 

A(I + h n A$ 1 S) = h n &iS . (4.12) 

Hence, (|4.5|) follows provided that I + /inA^iS 1 is invertible. □ 



Therefore any discrete gradient scheme has a locally exact modification 
of the form f 14 . 3 1) provided that the right-hand side of (J4.5P exists. This 
modification is unique if F' is invertible. By (j2.4]l ipi(h n F') is invertible 
for sufficiently small /i n and <^i(0) = J. Hence, (/ + h n A§iS)~ exists for 
sufficiently small h n and 

lim — A = S . (4.13) 

h„-t0 h n 

It turns out that A given by ( 14.51) is skew- symmetric, A T = —A. We point 
out that the skew symmetry of A is not assumed and is not obvious. In order 
to prove this fact we need the following lemma. 

Lemma 4.3. Suppose g(z) is an even analytic function (g(—z) = g(z)) 
and M is a matrix admitting the factorization M = QT , where T T = T , 
Q T = —Q and Q is invertible. Then 

(g(M)f = Q- 1 g(M)Q . (4.14) 



Proof: We represent gas a series 

oo 

g(M) = J2*k(QT) 2k . (4.15) 

fc=0 

Using assumed properties of Q, T we have (QT) = —TQ. Therefore 

oo oo 

(g(M)f = J2a k (TQ) 2k = J2^Q- 1 (QT) 2h Q . (4.16) 

fc=0 fc=0 

Hence, (|4~T4"1) follows. □ 



Theorem 4.4. Matrix A given by ( [^.5] ) is skew-symmet 



ric. 



Proof: We assume that h n is sufficiently small so that 3>i and A are both invertible. 
Then, (|4TT2"j) implies 

h n k- 1 = S- 1 ($^ + KSA) . (4.17) 

We define R = A — B. Then, taking into account SA + SB = F', we get 

SA = -F' + -SR , SB = -F' - -SR . (4.18) 

2 2' 22 v ' 

Hence, 

hnA- 1 = hx n R + S- x g{F') , 5 (F') = ^ + hi n F' . (4.19) 

Note that g(z) is an analytic function on the disc h n \z\ < 2ir, see (|2.4|) . Using 
Theorem 13.31 we easily show skew-symmetry of R: 

R T = A T - B T = B - A = -R . (4.20) 

Moreover, expressing ipi(h n F') is terms of e hnF (compare (|2.5p ). we have 

g(F') = (tpiihnF'))- 1 + ~h n F> = ±h n F' coth {^KF^ , (4.21) 

where zcoth(z) = zcosh(z)/sinh(z) is analytic for \z\ < ir. Hence g(-F') = g(F'). 
Moreover, F' = SH yy , where S T = —S and H^y = H yy . Therefore, by Lemma I3~3l 



{g(F')) T = S- l g(F')S . (4.22) 
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Then, using S 1 = S T and S 2 = —I, we obtain 

(S^g^f = (g(F')) T S = -S'^F') . (4.23) 

Taking into account (|4.20p and (|4.23p . we have from (J4.19P that A -1 is skew sym- 
metric. Hence A T = — A. □ 



Therefore A from ( 14. 5 p satisfies (14. 2p . which means that the locally ex- 
act modification defined by (14. 5 p is also energy preserving (for any discrete 
gradient VH). 



5 Special cases 

Locally exact modifications have simplest form in the case of symmetric 
discrete gradients. Then A = B and R = 0. In this case ( I4.19P and ( 14.2 ip 
yield: 

A = Manhc ( -h n F'\ S, (5.1) 

where tanhc(2:) := z _1 tanhz for z ^ and tanhc(O) := 1. The function 
tanhc(z) is analytic on the disc {z G C : \z\ < 7r/2}, hence A is well defined for 
sufficiently small h n , compare (12.61) . Note that formula ( 15. IP is independent 
of the discrete gradient for all symmetric discrete gradients. 

Specializing results of the previous section to separable Hamiltonians of 
the form 

H{x,p) = T(p) + V(x) (5.2) 

and to any symmetric discrete gradient V s , we get the following useful theo- 
rem, compare [12], Proposition 6.12. We use notation tanc(z) := z^ 1 tan(z) 
for z 7^ and tanc(0) := 1 (functions tanc(;z) and tanhc(^) are analytic in a 
neighbourhood of z — 0). 

Theorem 5.1. Numerical scheme 

8' 1 (x n+1 -x n ) = V s T(p n ,p n+ i) 

(5.3) 

( S l) * (Pn+1 -Pn) = ~ ^ s V(x n , X n+1 ) 
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preserves exactly the energy integral (for any h n - dependent invertible matrix 
S n such that S n — y h n I for h n — y 0) and is locally exact for 



n n \L n 



5 n = Mane—— , tt n = T„(p n )V m (x n ), (5.4) 

(S n given by \5.$ is invertible for sufficiently small h n , at least). 



Proof: The system (|5.3[) is equivalent to (|4.3|) if we take 

Hence, by Lemma 14. 11 the numerical scheme (|5.3p is energy preserving. 

Since V-ff is a symmetric discrete gradient if A is given by (|5.ip . then the 
scheme defined by (|5.3j) is locally exact. We have F = (T p , — V X ) T and 

F '=(-lt)> w— (^- ^J- <5 - 6) 

where all quantities are evaluated at (x n ,p n ). Denoting Q,"^ = TppV X x we have: 

? A2 f W 



<^-(o"(^J' (5 ' 7) 

Comparing Taylor series of tanhc(z) and tanc(z) we see that tanhc(A) = tanc(S) 
for any matrices satisfying B 2 = —A 2 . Therefore, formula (J5.7P implies 

**££-("*^ tmc (^))- <*■»> 

Therefore, if we define S n by (|5.4p . then the formula (|5.5p yields (|5.ip . which means 
that the scheme given by (j5.3f) is locally exact. □ 



In the case when x and p are scalars (x and p), then (-F') 2 is proportional 
to the unit matrix. Therefore, an energy-preserving locally exact scheme of 
the form 

X n +1 ~ %n &H p n+l - p n AH 



5 n Ap ' 5 n Ax 



(5.9) 



exists for any Hamiltonian and any symmetric discrete gradient, compare 
[TT] . It is enough to take 

h n u n 



8 n = /i n tanc— — , u n = JH XX H PP - H 2 p , (5.10) 
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where u n is evaluated at x n ,p n . 

We point out that the formula f)5. 10|) implies some limitations on h n in 
the case u n G R. Certainly we have to require h n u n ^ n + 2irM (M G N), or 
even h n u n < n. The last inequality is not very restrictive because it means 
that h n < |T n , where T n = 2Tr/u n is a corresponding period, see [TO]. 



The case H(x,p) = \p 2 + V(x) was considered in previous papers 
PHI [JTJ [TB] , where one can find results of many numerical experiments. As- 
suming x n = x n , p n = p n we obtain a scheme of 3rd order (GR-LEX), while 
\ [x n + x n+ i), p n = h (p n + Pn+i) yields a scheme of 4th order (GR- 



X i 



SLEX). In the case of one degree of freedom the discrete gradient scheme 
without modifications is of second order. Numerical experiments have shown 
that the accurcy of GR-LEX and GR-SLEX is higher by several orders of 
magnitude when compared with the standard discrete gradient method (while 
the computational cost is higher at most by several times, usually much less) 
[TO], ITT]. We point out, however, that in the multidimensional case the order 
of locally exact modifications is usually not greater than 2, and can be higher 
only for exceptional discrete gradients. 

In this Letter we present few results of numerical experiments for a two- 
dimensional anharmonic oscillator 

H(x,p)= 1 -\p\ 2 + V(\x\), y(|*|) = ±| a f--L|»|*. (5.11) 

see Figs. [TJ [2J OandSJ We consider circular orbits (of radius R), when the 
exact solution can be easily found (in this potential the radius of a circular 
orbit have to be smaller than 5). We compare the coordinate increment 
discrete gradient (J3.8P (denoted by GR), its symmetric modification (13. 9p 
(denoted by GR-SYM), and locally exact modifications ( 15. 3ft : GR-LEX (x = 
x n ) and GR-SLEX (x = \{x n +x n+ i)). In this case GR-SYM, GR-LEX and 
GR-SLEX are of the second order, while GR is only of the first order. 

Locally exact modifications are more expensive. The cost was estimated 
by the number of function evaluations. The average number of function 
evaluations per step depends on h and on the method. For instance, in the 
case of GR we have 85 evaluation per step for h = 0.05 and 160 evaluations 
per step for h = 0.5. In the case of GR-SLEX we have 262 and 341 evaluations 
per step, respectively. In numerical experiments we use different time steps 
in order to get the same computational costs. 

Fig. [TJ shows that for R = 0.1 locally exact modifications are more accu- 
rate than standard discrete gradient schemes by about 3 orders of magnitude. 
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We had to divide this figure into two parts (with different scales on axes). 
Figs. H] and [3] concern the case R = 1. Locally exact schemes are more ac- 
curate by several times. Note that both figures are quite similar although 
time steps are much smaller at Fig. |2j The efficiency of locally exact modi- 
fications decreases with increasing R. For R = 3 all four considered schemes 
have similar accuracy but surprisingly GR is the best, see Fig. HJ In the 
case of larger R our modifications are less accurate, especially when com- 
pared with GR-SYM. Therefore, we can conclude that locally exact schemes 
are very accurate in a neighbourhood (not very small, in fact) of the stable 
equilibrium. 

6 Concluding remarks 

We presented a construction of new numerical schemes based on the notion 
of local exactness. This notion has been known (under different names) for 
almost fifty years. The original application was confined to exact discretiza- 
tion of linearized equations [T3]. Our approach has two new features. First, 
we modify a given numerical scheme in a locally exact way. Second, we 
try to preserve geometric properties of the original numerical scheme. This 
task is not trivial. In this paper we present one successful application: lo- 
cally exact modifcations of discrete gradient methods for canonical Hamilton 
equations. We constructed a locally exact integrator (a modification of the 
discrete gradient scheme) which preserves exactly the energy integral. 

In the case of one degree of freedom the proposed modification, although 
more expensive (by only a few percent), turns out to be more accurate by as 
much as 8 orders of magnitude (in the case of small oscillations around the 
stable equilibrium) in comparison to the standard discrete gradient scheme 
[T0l[lT]. In multidimensional cases the relative cost of our algorithm is higher, 
but still we hope that our method will be of advantage. In general, the 
accuracy of locally exact algorithms is very high (much higher than their 
order suggests) in the neighbourhood of stable equilibria. Modifications 
proposed in this Letter contain exponentials of variable matrices. Similar 
time-consuming evaluations are characteristic for all exponential integrators. 
Effective methods of computing matrix exponentials, recently developed in 
this context [TH [29] , may decrease the computational cost of our algorithms. 

Another advantage of the proposed approach is the natural possibility of 
using a variable time step. Unlike symplectic methods (which work mostly 
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for constant time step) discrete gradient methods admit conservative modi- 
fications with variable time step. Therefore, one may easily implement any 
variable step method in order to obtain further improvement. 
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Figure 1: Error of numerical solutions for a circular orbit (R = 0.1, exact 
period T as 6.28) in the potential V(r) = 0.5r 2 — O.Olr 4 as a function of t. 
GR: dark line (h = 0.5), GR-SYM: black line (h = 0.625), GR-LEX: gray 
line (h = 0.766), GR-SLEX: light gray line (h = 1.063). 
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Figure 2: Error of numerical solutions for a circular orbit (R = 1.0, exact 
period T ps 6.41) in the potential V(r) = 0.5r 2 — O.Olr 4 as a function of t. 
GR: dark line (/i = 0.05), GR-SYM: black line (h = 0.067), GR-LEX: gray 
line (h = 0.094), GR-SLEX: light gray line (h = 0.154). 
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Figure 3: Error of numerical solutions (at t = 641) for a circular orbit (R = 
1.0, exact period T ss 6.41) in the potential V(r) = 0.5r 2 — O.Olr 4 as a 
function of t. GR: dark line (h = 0.5), GR-SYM: black line (h = 0.625), 
GR-LEX: gray line (h = 0.766), GR-SLEX: light gray line (h = 1.063). 
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Figure 4: Error of numerical solutions for a circular orbit (R = 3.0, exact 
period T w 7.85) in the potential V(r) = 0.5r 2 — O.Olr 4 as a function of t. 
GR: dark line (/i = 0.5), GR-SYM: black line (h = 0.627), GR-LEX: gray 
line (h = 0.768), GR-SLEX: light gray line (h = 1.066). 




20 



